Fix 1 (primary — caused the crash of test_mode_coeffs.py): Remove the geom_fix_object_ptr(&o) call from geom_get_bounding_box(). For Prism objects, this call triggered reinit_prism(), which frees and re-mallocs all internal prism arrays. When multiple OMP threads call this concurrently on the same Prism, it causes double-free/use-after-free, corrupting heap metadata and producing the free(): unaligned chunk detected in tcache 2 error.

This changes the semantics. But it's probably a good idea anyway since fix_object is expensive and should really only be done once at the beginning of all geometry calculations, and in fact we are already doing this.. (In Meep, it is called in the geom_epsilon constructor, and in MPB and pympb it is called in init_epsilon.

Change intersect_line_segment_with_prism() to use a stack-allocated variable-length array (VLA) (double slist[num_vertices + 2]) instead of the shared prsm->workspace.items buffer.

This is potentially dangerous for a large prism because stack space is limited. In particular, it could crash for prisms more than about 100,000 elements.

For the prism list, in practice the slist normally doesn't need to be equal to the number of vertices or anything close, since it is very unlikely you will have so many intersections. What we could do to make it both thread-safe and memory-safe would be something like:

1. Stack allocate a moderate number of intersections, e.g. double slist_stack[1024]; double *slist = slist_stack;, that will easily fit in the stack.
2. In intersect_line_with_prism, add an extra argument int slist_len for the length of slist (passing 1024)
3. Change the line from slist[num_intersections++] = s; to something like if (num_intersections++ < slist_len) slist[num_intersections-1] = s. Then, at the end of that loop, add a check if (num_intersections > slist_len) return num_intersections; to return immediately
4. In the caller intersect_line_segment_with_prism, check if (num_intersections > 1024) { slist = (double *) malloc(num_intersections * sizeof(double)); num_intersections = intersect_line_with_prism(prsm, pc, dc, slist, num_intersections); } to re-do the intersect_line_with_prism calculation with a heap-allocated array, as a slow fallback path. At the end of the function, also add if (num_intersections > 1024) free(slist); to deallocated it in this case.
5. Get rid of prsm->workspace, which is no longer used.

I implemented my proposal and pushed a commit to this PR.

The latest changes (c617512) are producing failures for the Meep unit test test_mode_coeffs.py when OMP_NUM_THREADS is set to 2 or 4:

OMP_NUM_THREADS=2:

malloc(): unaligned tcache chunk detected

OMP_NUM_THREADS=4:

degenerate plane in project_point_into_plane